OVERVIEW OF MAIN STEPS IN RNA-Seq DATA ANALYSIS

(da inserire tabella)

RAW READS MANIPULATION

Accendere il computer e aprire il terminale (ctrl + alt + T)

Spostarsi nella cartella Scaricati, create una nuova cartella chiamandola tc

cd Scaricati/
mkdir tc

Entrare nella cartella appena creata e valutarne il contenuto della cartella digitare il seguente comando

cd tc/
ls -lrth

La cartella dovrebbe essere vuota. Come prima cosa dobbiamo recuperare i file grezzi (FASTQ) del nostro esperimento. Per recuperare i file ci connettiamo al nostro server di laboratorio e con il protocollo di trasferimento ssh copiamo in locale (la cartella dove siamo in questo momento) i 2 file corrispondenti alle reads forward e alle reads reverse. Lavoreremo con un solo campione, ma come sapete l’esperimento consiste di 3 repliche biologiche per i non trattati e 3 repliche biologiche per quelli trattati. I comandi che qui eseguiremo possono essere ripetuti su più file contemporaneamente tramite appositi comandi come i cicli for

# prima scarichiamo il primo file (attenzione al punto in fondo al comando, indica current directory)
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_1_subset.fastq.gz .
# inserire la password quando ce la chiede ()
# controlliamo che sia stato scaricato
ls -lrth
# ora scarichiamo il secondo file
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/Sample_Untreated_1_2_subset.fastq.gz .
# inserire nuovamente la password
# controlliamo
ls -lrth

Come potete osservare la dimensione di questi file è molto bassa rispetto alla dimensione tipica di un file che risulta da un esperimento di sequenziamento come lo abbiamo descritto. Questo perché i file originali sono stati copiati e ridotti in dimensioni con lo scopo di riuscire a manipolarli e lavorarci nei tempi della lezione, principalmente per il fatto che le capacità dei computer che stiamo usando potrebbero essere limitate e quindi non in grado di eleaborare questi file in tempi brevi

Iniziamo con esplorare questi file

# "apriamo" il file
less -S Sample_Untreated_1_1_subset.fastq.gz
# contiamo il numero di reads
zcat Sample_Untreated_1_1_subset.fastq.gz | echo $((`wc -l`/4))
zcat Sample_Untreated_1_2_subset.fastq.gz | echo $((`wc -l`/4))
# il numero di reads F e R coincide?

QUALITY CONTROL

Per effettuare il controllo di qualità sulle nostre reads utilizzeremo il programma FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Questo software (scritto in Java) prende come input un numero a piacimento di reads e calcola varie statistiche che possono dare un’idea della qualità dei nostri dati

# creiamo una nuova cartella dove scaricheremo il programma e salveremo i risultati
mkdir fastqc
cd fastqc/
# scarichiamo il programma con il comando wget
wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
# scompattiamo la cartella
unzip fastqc_v0.11.7.zip
# rimuoviamo cartella compressa
rm fastqc_v0.11.7.zip
# vediamo cosa c'è dentro la cartella senza spostarci dentro (percorso relativo alla nostra posizione)
ls -lrth FastQC/
# il file che ci interessa, cioè il programma che andremo a richiamare è "fastqc"
# dobbiamo prima renderlo eseguibile
chmod 755 FastQC/fastqc
# controlliamo che funzioni
./FastQC/fastcq -h

A questo punto possiamo lanciare il programma dando come input i due file con le reads F e R. Il comando è molto semplice, basta specificare i nomi dei file che vogliamo controllare, e la cartella dove vogliamo che vengano salvati i report coi risultati (la cartella deve già esistere quindi prima la creiamo)

# creiamo cartella per salvare i risultati
mkdir results
# ci spostiamo una cartella su, dove stanno i nostri file
cd ..
# eseguiamo il programma (dobbiamo specificare il percorso dell'eseguibile rispetto a dove siamo) e specifichiamo inoltre il percorso della cartella di output
./fastqc/FastQC/fastqc Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz -o fastqc/results/

Per ogni file di input vengono creati 2 file di output, un file html e una cartella compressa. Andiamo ad esplorare il report in formato html aprendolo con Firefox

FastQC report

FastQC report

A questo punto commentiamo insieme i risultati del controllo qualità confrontandoci con qualità di riferimento sul sito di FastQC e con la spiegazione dei moduli di analisi (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/)

Dopo aver deciso insieme quali parametri utilizzare per correggere la qualità delle nostre reads ed aver identificato gli adattatori da rimuovere, andiamo ad effettuare la rimozione degli adattatori e il quality trimming

Adapters removal and quality trimming

Per effettuare queste due operazioni useremo un unico programma che si chiama Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). Il suo utilizzo è molto simile a quello di FastQC (anche Trimmomatic è scritto in Java) e necessita come input le sequenze grezze. Per ogni file di input (nel nostro caso 2 file) verranno creati 2 file di output, uno indicato come paired, l’altro come unpaired; è necessario indicare il nome di ciascuno dei file di output. Le opzioni di Trimmomatic sono:

# ci spostiamo nella cartella tc, creiamo nuova cartella e scarichiamo il programma
# usiamo percorso assoluto da qualsiasi punto del computer per spostarci dentro una qualsiasi cartella
cd /home/usr/Scaricati/tc # usr è specifico per il vostro utente (il mio è dlfptr70)
mkdir Trimmomatic
wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.36.zip -O Trimmomatic/Trimmomatic-0.36.zip # con l'opzione -O specifichiamo dove vogliamo scaricare il file e con che nome
# decomprimiamo senza cambiare cartella e rimuoviamo l'archivio
unzip Trimmomatic/Trimmomatic-0.36.zip -d Trimmomatic/ && rm Trimmomatic/Trimmomatic-0.36.zip
# controlliamo il funzionamento del programma
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar -h # -h è l'opzione per help

A questo punto possiamo lanciare Trimmomatic sui nostri file; anche se dal controllo qualità non abbiamo identificato sequenze di adattatori come contaminanti, un classico approccio potrebbe essere quello di dare come input di file degli adattatori il file fornito dal programma che quindi rimuoverà tutte le sequenze corrispondenti a quelle che trova

# creiamo file unendo in uno tutti gli adattatori
cat Trimmomatic/Trimmomatic-0.36/adapters/*.fa > Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa
# cat è un comando che unisce file "row-wise"; l'asterisco (*.fa) indica di selezionare tutti i file che hanno estensione .fa
java -jar Trimmomatic/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 10 -trimlog \
Trimmomatic/log_trimmomatic \
Sample_Untreated_1_1_subset.fastq.gz Sample_Untreated_1_2_subset.fastq.gz \
Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_1_subset_trim_unpaired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz Trimmomatic/Sample_Untreated_1_2_subset_trim_unpaired.fastq.gz \
ILLUMINACLIP:Trimmomatic/Trimmomatic-0.36/adapters/ALL_ADAPTERS.fa:2:30:10 \
LEADING:16 TRAILING:16 SLIDINGWINDOW:4:25 MINLEN:36 &> Trimmomatic/log_trimmomatic_screen
# con il simbolo &> andiamo a salvare su un file anche l'output che normalmente sarebbe stampato sullo schermo

Dovrebbe terminare dopo qualche secondo senza aver stampato a schermo alcuna informazione. Andiamo a vedere se ha creato i file ed esaminiamo i file di log (log_trimmomatic e log_trimmomatic_screen); il primo contiene le seguenti info per ciascuna read:

Il secondo contiene le statistiche generali

# guardiamo prima le statistiche generali
less -S Trimmomatic/log_trimmomatic_screen
# poi l'altro file
less -S Trimmomatic/log_trimmomatic

A questo punto ripetiamo il controllo di qualità, questa volta sulle reads trimmate, solo sui file paired

# creiamo cartella dove salveremo nuovo controllo qualità
mkdir fastqc/results_after_trimming
# lanciamo FastQC sulle sequenze trimmate
./fastqc/FastQC/fastqc Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz -o fastqc/results_after_trimming/

Andiamo ad esplorare i report e vedere se le statistiche sulla qualità sono cambiate

MAPPING (ALLINEAMENTO)

Ora i nostri file sono pronti per essere allineati sul genoma di riferimento (Arabidopsis). Per quasi tutti i programmi che effettuano allineamento il primo passaggio consiste della creazione dell’indice del genoma. Perché viene creato l’indice del genoma: citazione presa da qualcuno ma non ricordo chi (modificata) “Se entro in una biblioteca con la frase “La genetica è il ramo della biologia che si occupa del materiale ereditario” e voglio trovare il libro in cui è scritta, farò molta fatica a trovarlo, a meno che all’ingresso della biblioteca non si trovi un catalogo che mi indirizza almeno verso la sezione “Scienza””. L’indice quindi consiste in una serie di file che riassumono le principali caratteristiche del genoma di un organismo. Per effettuare l’allineamento useremo un programma chiamato STAR (https://github.com/alexdobin/STAR). Andiamo a scaricarlo insieme ai file aggiornati del genoma e del trascrittoma di Arabidopsis

# siamo nella cartella tc
mkdir STAR
# scarichiamo e salviamo nella cartella appena creata
wget https://github.com/alexdobin/STAR/archive/2.6.0c.tar.gz -O STAR/2.6.0c.tar.gz
# decomprimiamo e rimuoviamo archivio
tar -xf STAR/2.6.0c.tar.gz -C STAR/ && rm 2.6.0c.tar.gz
# controlliamo il funzionamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR -h

Per creare l’indice del genoma, abbiamo bisogno del genoma. Andiamo a scaricarlo dal sito di Ensembl Plants

# genoma
mkdir genome
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz -O genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
# trascrittoma
wget ftp://ftp.ensemblgenomes.org/pub/release-39/plants/gtf/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.39.gtf.gz -O Arabidopsis_thaliana.TAIR10.39.gtf.gz
# i file devono essere decompressi
gunzip genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip genome/Arabidopsis_thaliana.TAIR10.39.gtf.gz

A questo punto possiamo creare l’indice con il comando di STAR --runMode genomeGenerate. Vediamo nel dettaglio le altre opzioni:

Lanciamo il comando con le opzioni e i percorsi adattati al nostro caso

# creiamo cartella per l'indice all'interno della cartella del genoma
mkdir -p genome/index
# lanciamo la creazione dell'indice
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--runMode genomeGenerate \
--genomeDir genome/index \
--genomeFastaFiles genome/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa \
--sjdbGTFfile genome/Arabidopsis_thaliana.TAIR10.39.gtf \
--sjdbOverhang 100

Con 10 core dovrebbe impiegare circa 3 minuti. Andiamo a vedere se ha creato i file nella cartella index

# spostiamo il file Log.out dove dovrebbe stare
mv Log.out genome/index/
# guardiamo i file appena creati
ls -lrht/genome/index
# raccolgono diverse statistiche sul genoma di Arabidopsis, comprese info sul db di splicejunctions

Possiamo ora andare ad effettuare l’allineamento. Ovviamente per l’allineamento useremo le reads trimmate e pulite (le paired). Vediamo nel dettaglio i comandi:

Attenzione che la maggior parte dei programmi di allineamento sono creati con impostazioni di default basate sul genoma umano –> aggiustare i parametri sulla base del proprio organismo! Andiamo a lanciare il comando

# creiamo cartella dove salveremo file di allineamento
mkdir -p STAR/aligned
# lanciamo l'allineamento
./STAR/STAR-2.6.0c/bin/Linux_x86_64_static/STAR \
--runThreadN 10 \
--genomeDir genome/index/ \
--readFilesIn Trimmomatic/Sample_Untreated_1_1_subset_trim_paired.fastq.gz \
Trimmomatic/Sample_Untreated_1_2_subset_trim_paired.fastq.gz \
--readFilesCommand zcat \
--outFileNamePrefix STAR/aligned/SU_1_ \
--alignIntronMax 2000 \
--outSAMtype BAM SortedByCoordinate

Se è tutto a posto l’allineamento dovrebbe impiegare circa 10 secondi. Come vedete la lunghezza massima media per gli introni di Arabidopsis è di circa 2000 bp, e gli ho detto che voglio che il file di output sia in formato .bam (formato allineamento binario). Il file .bam è in un formato che non può essere aperto con il terminale, ma che andremo poi a visualizzare con un apposito software. Per completezza vi inserisco un immagine dello stesso file di allineamento ma in un formato visualizzabile e consultabile sul terminale. Il formato di questo file è .sam (Sequence Alignment Format), formato non molto utilizzato perché non essendo compresso occupa molto spazio

SAM file

SAM file

Andiamo ora a vedere i nostri file creati

# il file che ci interessa è chiamato SU_1_Aligned.sortedByCoord.out.bam
ls -rht STAR/aligned/
# leggiamo qualche statistica in questo altro file
less -S STAR/aligned/SU_1_Log.final.out

Quanto vi risulta la percentuale di Uniquely mapped reads?

ALIGNMENT VISUALIZATION

Dopo aver effettuato l’allineamento può essere utile effettuare un controllo visivo sull reads allineate. A tal fine si utilizza spesso un software chiamato IGV (http://software.broadinstitute.org/software/igv/) che permette tra le altre cose di visualizzare i file dell’allineamento dopo aver caricato un genoma di riferimento

Integrative Genomics Viewer

Integrative Genomics Viewer

Andiamo a scaricare il programma dal sito http://software.broadinstitute.org/software/igv/download cliccando in basso su “Launch with 2 GB”. Attendiamo lo scaricamento e apriamo il programma con doppio clic. Una volta aperto in alto a sinistra selezionare il genoma di A. thaliana TAIR10. A questo punto possiamo caricare i risultati del nostro allineamento, non prima però di aver creato un indice per il file .bam. Un altro indice direte voi, sì questa volta però non eseguiremo nessun passaggio ma vi farò prelevare il file dell’indice già preparato dal nostro server. Per completezza vi mostro il comando utilizzato per creare l’indice con Samtools

samtools index -b SU_1_Aligned.sortedByCoord.out.bam

Su questi computer non possiamo installare il programma che si utilizza, tra le altre cose, anche per la creazione degli indici per i file di allineamento. Il programma si chiama Samtools (http://www.htslib.org/) ed è molto utilizzato in bioinformatica perché svolge numerose funzioni. Andiamo a copiare nella cartella dove sta il file .bam appena creato da noi, il file dell’indice, che permette che il file .bam venga visualizzato correttamente in IGV

scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/
# il file che vedete deve avere lo stesso nome del file di allineamento, ma ha un'estensione diversa, .bai

A questo punto potete tornare su IGV e cliccare su “File -> Load from file” e selezionare dalla cartella aligned il file SU_1_Aligned.sortedByCoord.out.bam. Selezioniamo Chr1 e zoomiamo sino a che la rotella inizia a girare, a quel punto possiamo visualizzare le nostre reads allineate ai geni del genoma di riferimento

Non sarà così facile trovare un po’ di reads!

Non sarà così facile trovare un po’ di reads!

Nel caso non funzionasse, copiare dal server alla vostra cartella, sia il file .bam che il .bai (potrebbe essere che il bai debba essere creato dallo stesso bam per funzionare la visualizzazione)

# bam
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam STAR/aligned/
# bai
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/SU_1_Aligned.sortedByCoord.out.bam.bai STAR/aligned/

La lunghezza degli introni è stata rispettata?

READS COUNTING

Una volta controllato che l’allineamento è andato a buon fine e siamo soddisfatti di come il programma di mapping ha lavorato, possiamo spostarci nel prossimo step verso l’analisi dei geni, ovvero la conta delle reads. Per conta delle reads si intende il contare quante reads sono state assegnate (hanno “mappato”) a un determinato gene (o trascritto, o esone, le cosiddette genomic features) in seguito al mapping. Per questo passaggio utilizzeremo il software featureCounts (http://subread.sourceforge.net/) del pacchetto Subread. Scarichiamo e testiamo il programma

# creiamo cartella
mkdir featureCounts
# scarichiamo dal server
scp trascrittomica@157.27.231.104:/home/trascrittomica/Scaricati/featureCounts featureCounts/
# testiamo funzionamento
./featureCounts/featureCounts -v

Vediamo le opzioni e i parametri di featureCounts:

A questo punto andiamo a lanciare il programma, ricordandoci che nel nostro caso vogliamo contare le feature del file di allineamento .bam

# creiamo cartella per output
mkdir -p featureCounts/output
# lanciamo la conta
./featureCounts/featureCounts \
-a genome/Arabidopsis_thaliana.TAIR10.39.gtf \
-o featureCounts/output/SU_1_counts.txt \
-p \
-t exon \
-g gene_id \
-T 10 \
STAR/aligned/SU_1_Aligned.sortedByCoord.out.bam

Dovrebbe terminare in qualche secondo. Andiamo ad esplorare i file di output

ls -lrht featureCounts/output/
less -S featureCounts/output/SU_1_counts.txt.summary
less -S featureCounts/output/SU_1_counts.txt

Le colonne che ci interessano sono la prima e l’ultima (7)